{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "using SymPy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "E (generic function with 2 methods)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Create an n x n (default 3), i,j elimination matrix\n",
    "function E(i,j,n=3)\n",
    "    E = sympy\"eye\"(n)           # start with the identity\n",
    "    E[i,j] = - symbols(\"m$i$j\")  # insert the negative multiplier\n",
    "    E\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{bmatrix}A_{11}&A_{12}&A_{13}\\\\A_{21}&A_{22}&A_{23}\\\\A_{31}&A_{32}&A_{33}\\end{bmatrix}"
      ],
      "text/plain": [
       "3×3 Array{SymPy.Sym,2}:\n",
       " A11  A12  A13\n",
       " A21  A22  A23\n",
       " A31  A32  A33"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Create a symbolic 3x3 A matrix\n",
    "  A = [symbols(\"A$i$j\") for i=1:3, j=1:3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{bmatrix}1&0&0\\\\- m_{21}&1&0\\\\0&0&1\\end{bmatrix}"
      ],
      "text/plain": [
       "3×3 Array{SymPy.Sym,2}:\n",
       "    1  0  0\n",
       " -m21  1  0\n",
       "    0  0  1"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E(2,1) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div id=\"interact-js-shim\">\n",
       "    <script charset=\"utf-8\">\n",
       "(function (IPython, $, _, MathJax) {\n",
       "    $.event.special.destroyed = {\n",
       "\tremove: function(o) {\n",
       "\t    if (o.handler) {\n",
       "\t\to.handler.apply(this, arguments)\n",
       "\t    }\n",
       "\t}\n",
       "    }\n",
       "\n",
       "    var OutputArea = IPython.version >= \"4.0.0\" ? require(\"notebook/js/outputarea\").OutputArea : IPython.OutputArea;\n",
       "\n",
       "    var redrawValue = function (container, type, val) {\n",
       "\tvar selector = $(\"<div/>\");\n",
       "\tvar oa = new OutputArea(_.extend(selector, {\n",
       "\t    selector: selector,\n",
       "\t    prompt_area: true,\n",
       "\t    events: IPython.events,\n",
       "\t    keyboard_manager: IPython.keyboard_manager\n",
       "\t})); // Hack to work with IPython 2.1.0\n",
       "\n",
       "\tswitch (type) {\n",
       "\tcase \"image/png\":\n",
       "            var _src = 'data:' + type + ';base64,' + val;\n",
       "\t    $(container).find(\"img\").attr('src', _src);\n",
       "\t    break;\n",
       "\tcase \"text/latex\":\n",
       "\t\tif (MathJax){\n",
       "\t\t\tvar math = MathJax.Hub.getAllJax(container)[0];\n",
       "\t\t\tMathJax.Hub.Queue([\"Text\", math, val.replace(/^\\${1,2}|\\${1,2}$/g, '')]);\n",
       "\t\t\tbreak;\n",
       "\t\t}\n",
       "\tdefault:\n",
       "\t    var toinsert = OutputArea.append_map[type].apply(\n",
       "\t\toa, [val, {}, selector]\n",
       "\t    );\n",
       "\t    $(container).empty().append(toinsert.contents());\n",
       "\t    selector.remove();\n",
       "\t}\n",
       "    }\n",
       "\n",
       "\n",
       "    $(document).ready(function() {\n",
       "\tfunction initComm(evt, data) {\n",
       "\t    var comm_manager = data.kernel.comm_manager;\n",
       "        //_.extend(comm_manager.targets, require(\"widgets/js/widget\"))\n",
       "\t    comm_manager.register_target(\"Signal\", function (comm) {\n",
       "            comm.on_msg(function (msg) {\n",
       "                var val = msg.content.data.value;\n",
       "                $(\".signal-\" + comm.comm_id).each(function() {\n",
       "                var type = $(this).data(\"type\");\n",
       "                if (typeof(val[type]) !== \"undefined\" && val[type] !== null) {\n",
       "                    redrawValue(this, type, val[type], type);\n",
       "                }\n",
       "                });\n",
       "                delete val;\n",
       "                delete msg.content.data.value;\n",
       "            });\n",
       "\t    });\n",
       "\n",
       "\t    // coordingate with Comm and redraw Signals\n",
       "\t    // XXX: Test using Reactive here to improve performance\n",
       "\t    $([IPython.events]).on(\n",
       "\t\t'output_appended.OutputArea', function (event, type, value, md, toinsert) {\n",
       "\t\t    if (md && md.reactive) {\n",
       "                // console.log(md.comm_id);\n",
       "                toinsert.addClass(\"signal-\" + md.comm_id);\n",
       "                toinsert.data(\"type\", type);\n",
       "                // Signal back indicating the mimetype required\n",
       "                var comm_manager = IPython.notebook.kernel.comm_manager;\n",
       "                var comm = comm_manager.comms[md.comm_id];\n",
       "                comm.then(function (c) {\n",
       "                    c.send({action: \"subscribe_mime\",\n",
       "                       mime: type});\n",
       "                    toinsert.bind(\"destroyed\", function() {\n",
       "                        c.send({action: \"unsubscribe_mime\",\n",
       "                               mime: type});\n",
       "                    });\n",
       "                })\n",
       "\t\t    }\n",
       "\t    });\n",
       "\t}\n",
       "\n",
       "\ttry {\n",
       "\t    // try to initialize right away. otherwise, wait on the status_started event.\n",
       "\t    initComm(undefined, IPython.notebook);\n",
       "\t} catch (e) {\n",
       "\t    $([IPython.events]).on('kernel_created.Kernel kernel_created.Session', initComm);\n",
       "\t}\n",
       "    });\n",
       "})(IPython, jQuery, _, MathJax);\n",
       "</script>\n",
       "    <script>\n",
       "        window.interactLoadedFlag = true\n",
       "       $(\"#interact-js-shim\").bind(\"destroyed\", function () {\n",
       "           if (window.interactLoadedFlag) {\n",
       "               console.warn(\"JavaScript required by Interact will be removed if you remove this cell or run using Interact more than once.\")\n",
       "           }\n",
       "       })\n",
       "       $([IPython.events]).on(\"kernel_starting.Kernel kernel_restarting.Kernel\", function () { window.interactLoadedFlag = false })\n",
       "   </script>\n",
       "</div>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "using Interact"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ca73636d-e2de-40e3-bf33-af0aefec0cb4",
       "version_major": 2,
       "version_minor": 0
      }
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [],
      "text/plain": [
       "Interact.Options{:SelectionSlider,Any}(1: \"input\" = 2 Any , \"i\", 2, \"2\", 2, Interact.OptionDict(DataStructures.OrderedDict{Any,Any}(\"1\"=>1,\"2\"=>2,\"3\"=>3), Dict{Any,Any}(Pair{Any,Any}(2, \"2\"),Pair{Any,Any}(3, \"3\"),Pair{Any,Any}(1, \"1\"))), Any[], Any[], true, \"horizontal\")"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c223c5da-a42c-4d7f-be1d-acea2be9152c",
       "version_major": 2,
       "version_minor": 0
      }
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [],
      "text/plain": [
       "Interact.Options{:SelectionSlider,Any}(3: \"input-2\" = 2 Any , \"j\", 2, \"2\", 2, Interact.OptionDict(DataStructures.OrderedDict{Any,Any}(\"1\"=>1,\"2\"=>2,\"3\"=>3), Dict{Any,Any}(Pair{Any,Any}(2, \"2\"),Pair{Any,Any}(3, \"3\"),Pair{Any,Any}(1, \"1\"))), Any[], Any[], true, \"horizontal\")"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "\\begin{bmatrix}1&0&0\\\\0&- m_{22}&0\\\\0&0&1\\end{bmatrix}"
      ],
      "text/plain": [
       "3×3 Array{SymPy.Sym,2}:\n",
       " 1     0  0\n",
       " 0  -m22  0\n",
       " 0     0  1"
      ]
     },
     "execution_count": 7,
     "metadata": {
      "comm_id": "e2514bf5-28a5-4eed-90ac-ba5ddfda6c73",
      "reactive": true
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@manipulate for i=1:3, j=1:3\n",
    "     E(i,j)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{bmatrix}1&0&0\\\\- m_{21}&1&0\\\\0&0&1\\end{bmatrix}"
      ],
      "text/plain": [
       "3×3 Array{SymPy.Sym,2}:\n",
       "    1  0  0\n",
       " -m21  1  0\n",
       "    0  0  1"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E(2,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{bmatrix}1&0&0\\\\m_{21}&1&0\\\\0&0&1\\end{bmatrix}"
      ],
      "text/plain": [
       "3×3 Array{SymPy.Sym,2}:\n",
       "   1  0  0\n",
       " m21  1  0\n",
       "   0  0  1"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "inv(E(2,1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "subtract m21 times the first row from the second row"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{bmatrix}A_{11}&A_{12}&A_{13}\\\\- A_{11} m_{21} + A_{21}&- A_{12} m_{21} + A_{22}&- A_{13} m_{21} + A_{23}\\\\A_{31}&A_{32}&A_{33}\\end{bmatrix}"
      ],
      "text/plain": [
       "3×3 Array{SymPy.Sym,2}:\n",
       "            A11             A12             A13\n",
       " -A11*m21 + A21  -A12*m21 + A22  -A13*m21 + A23\n",
       "            A31             A32             A33"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E(2,1) * A   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "do the row operation twice"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{bmatrix}A_{11}&A_{12}&A_{13}\\\\- 2 A_{11} m_{21} + A_{21}&- 2 A_{12} m_{21} + A_{22}&- 2 A_{13} m_{21} + A_{23}\\\\A_{31}&A_{32}&A_{33}\\end{bmatrix}"
      ],
      "text/plain": [
       "3×3 Array{SymPy.Sym,2}:\n",
       "              A11               A12               A13\n",
       " -2*A11*m21 + A21  -2*A12*m21 + A22  -2*A13*m21 + A23\n",
       "              A31               A32               A33"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E(2,1)^2 * A  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " subtract m32 times the second row from the third row"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{bmatrix}A_{11}&A_{12}&A_{13}\\\\A_{21}&A_{22}&A_{23}\\\\- A_{21} m_{32} + A_{31}&- A_{22} m_{32} + A_{32}&- A_{23} m_{32} + A_{33}\\end{bmatrix}"
      ],
      "text/plain": [
       "3×3 Array{SymPy.Sym,2}:\n",
       "            A11             A12             A13\n",
       "            A21             A22             A23\n",
       " -A21*m32 + A31  -A22*m32 + A32  -A23*m32 + A33"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E(3,2) * A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that after computing E(3,2) * A, the first row is untouched so we can apply E(2,1) without any interference from row 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{bmatrix}A_{11}&A_{12}&A_{13}\\\\- A_{11} m_{21} + A_{21}&- A_{12} m_{21} + A_{22}&- A_{13} m_{21} + A_{23}\\\\- A_{21} m_{32} + A_{31}&- A_{22} m_{32} + A_{32}&- A_{23} m_{32} + A_{33}\\end{bmatrix}"
      ],
      "text/plain": [
       "3×3 Array{SymPy.Sym,2}:\n",
       "            A11             A12             A13\n",
       " -A11*m21 + A21  -A12*m21 + A22  -A13*m21 + A23\n",
       " -A21*m32 + A31  -A22*m32 + A32  -A23*m32 + A33"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E(2,1) * E(3,2)  * A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, in the below, row 2 has changed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{bmatrix}A_{11}&A_{12}&A_{13}\\\\- A_{11} m_{21} + A_{21}&- A_{12} m_{21} + A_{22}&- A_{13} m_{21} + A_{23}\\\\A_{31}&A_{32}&A_{33}\\end{bmatrix}"
      ],
      "text/plain": [
       "3×3 Array{SymPy.Sym,2}:\n",
       "            A11             A12             A13\n",
       " -A11*m21 + A21  -A12*m21 + A22  -A13*m21 + A23\n",
       "            A31             A32             A33"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E(2,1) * A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "so if apply E(3,2) <br>\n",
    "(meaning:  subtract m32 times the second row from the third row) <br>\n",
    "it will happen with the UPDATED row 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{bmatrix}A_{11}&A_{12}&A_{13}\\\\- A_{11} m_{21} + A_{21}&- A_{12} m_{21} + A_{22}&- A_{13} m_{21} + A_{23}\\\\A_{11} m_{21} m_{32} - A_{21} m_{32} + A_{31}&A_{12} m_{21} m_{32} - A_{22} m_{32} + A_{32}&A_{13} m_{21} m_{32} - A_{23} m_{32} + A_{33}\\end{bmatrix}"
      ],
      "text/plain": [
       "3×3 Array{SymPy.Sym,2}:\n",
       "                         A11  …                          A13\n",
       "              -A11*m21 + A21                  -A13*m21 + A23\n",
       " A11*m21*m32 - A21*m32 + A31     A13*m21*m32 - A23*m32 + A33"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E(3,2) * E(2,1) * A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The row interpretation of matrix muliply correctly acounts for m32 times the second row and m21 x m32 times the first row -- the combined effect"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{bmatrix}1&0&0\\\\- m_{21}&1&0\\\\m_{21} m_{32}&- m_{32}&1\\end{bmatrix}"
      ],
      "text/plain": [
       "3×3 Array{SymPy.Sym,2}:\n",
       "       1     0  0\n",
       "    -m21     1  0\n",
       " m21*m32  -m32  1"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E(3,2) * E(2,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{bmatrix}1&0&0\\\\- m_{21}&1&0\\\\0&- m_{32}&1\\end{bmatrix}"
      ],
      "text/plain": [
       "3×3 Array{SymPy.Sym,2}:\n",
       "    1     0  0\n",
       " -m21     1  0\n",
       "    0  -m32  1"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E(2,1) * E(3,2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's move to 5x5 matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{bmatrix}1&0&0&0&0\\\\- m_{21}&1&0&0&0\\\\0&0&1&0&0\\\\0&0&0&1&0\\\\0&0&0&0&1\\end{bmatrix}"
      ],
      "text/plain": [
       "5×5 Array{SymPy.Sym,2}:\n",
       "    1  0  0  0  0\n",
       " -m21  1  0  0  0\n",
       "    0  0  1  0  0\n",
       "    0  0  0  1  0\n",
       "    0  0  0  0  1"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E(2,1,5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{bmatrix}1&0&0&0&0\\\\- m_{21}&1&0&0&0\\\\- m_{31}&0&1&0&0\\\\- m_{41}&0&0&1&0\\\\- m_{51}&0&0&0&1\\end{bmatrix}"
      ],
      "text/plain": [
       "5×5 Array{SymPy.Sym,2}:\n",
       "    1  0  0  0  0\n",
       " -m21  1  0  0  0\n",
       " -m31  0  1  0  0\n",
       " -m41  0  0  1  0\n",
       " -m51  0  0  0  1"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E1 = E(2,1,5) * E(3,1,5) * E(4,1,5)* E(5,1,5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{bmatrix}1&0&0&0&0\\\\- m_{21}&1&0&0&0\\\\- m_{31}&0&1&0&0\\\\- m_{41}&0&0&1&0\\\\- m_{51}&0&0&0&1\\end{bmatrix}"
      ],
      "text/plain": [
       "5×5 Array{SymPy.Sym,2}:\n",
       "    1  0  0  0  0\n",
       " -m21  1  0  0  0\n",
       " -m31  0  1  0  0\n",
       " -m41  0  0  1  0\n",
       " -m51  0  0  0  1"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E(3,1,5) * E(2,1,5) * E(5,1,5)* E(4,1,5)  # Why does the order not matter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{bmatrix}1&0&0&0&0\\\\- m_{21}&1&0&0&0\\\\- m_{31}&0&1&0&0\\\\- m_{41}&0&0&1&0\\\\- m_{51}&0&0&0&1\\end{bmatrix}"
      ],
      "text/plain": [
       "5×5 Array{SymPy.Sym,2}:\n",
       "    1  0  0  0  0\n",
       " -m21  1  0  0  0\n",
       " -m31  0  1  0  0\n",
       " -m41  0  0  1  0\n",
       " -m51  0  0  0  1"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E1 # subtracts multiples of the first row"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{bmatrix}1&0&0&0&0\\\\m_{21}&1&0&0&0\\\\m_{31}&0&1&0&0\\\\m_{41}&0&0&1&0\\\\m_{51}&0&0&0&1\\end{bmatrix}"
      ],
      "text/plain": [
       "5×5 Array{SymPy.Sym,2}:\n",
       "   1  0  0  0  0\n",
       " m21  1  0  0  0\n",
       " m31  0  1  0  0\n",
       " m41  0  0  1  0\n",
       " m51  0  0  0  1"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "inv(E1) # INVERSE means add back the same multiples of the first row"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{bmatrix}1&0&0&0&0\\\\0&1&0&0&0\\\\0&- m_{32}&1&0&0\\\\0&0&0&1&0\\\\0&0&0&0&1\\end{bmatrix}"
      ],
      "text/plain": [
       "5×5 Array{SymPy.Sym,2}:\n",
       " 1     0  0  0  0\n",
       " 0     1  0  0  0\n",
       " 0  -m32  1  0  0\n",
       " 0     0  0  1  0\n",
       " 0     0  0  0  1"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E(3,2,5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{bmatrix}1&0&0&0&0\\\\0&1&0&0&0\\\\0&m_{32}&1&0&0\\\\0&0&0&1&0\\\\0&0&0&0&1\\end{bmatrix}"
      ],
      "text/plain": [
       "5×5 Array{SymPy.Sym,2}:\n",
       " 1    0  0  0  0\n",
       " 0    1  0  0  0\n",
       " 0  m32  1  0  0\n",
       " 0    0  0  1  0\n",
       " 0    0  0  0  1"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "inv(E(3,2,5))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{bmatrix}1&0&0&0&0\\\\- m_{21}&1&0&0&0\\\\- m_{31}&- m_{32}&1&0&0\\\\- m_{41}&0&0&1&0\\\\- m_{51}&0&0&0&1\\end{bmatrix}"
      ],
      "text/plain": [
       "5×5 Array{SymPy.Sym,2}:\n",
       "    1     0  0  0  0\n",
       " -m21     1  0  0  0\n",
       " -m31  -m32  1  0  0\n",
       " -m41     0  0  1  0\n",
       " -m51     0  0  0  1"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E1 * E(3,2,5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{bmatrix}1&0&0&0&0\\\\m_{21}&1&0&0&0\\\\m_{31}&m_{32}&1&0&0\\\\m_{41}&0&0&1&0\\\\m_{51}&0&0&0&1\\end{bmatrix}"
      ],
      "text/plain": [
       "5×5 Array{SymPy.Sym,2}:\n",
       "   1    0  0  0  0\n",
       " m21    1  0  0  0\n",
       " m31  m32  1  0  0\n",
       " m41    0  0  1  0\n",
       " m51    0  0  0  1"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "inv(E1) * inv(E(3,2,5))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{bmatrix}1&0&0&0&0\\\\- m_{21}&1&0&0&0\\\\- m_{31}&- m_{32}&1&0&0\\\\- m_{41}&0&0&1&0\\\\- m_{51}&0&0&0&1\\end{bmatrix}"
      ],
      "text/plain": [
       "5×5 Array{SymPy.Sym,2}:\n",
       "    1     0  0  0  0\n",
       " -m21     1  0  0  0\n",
       " -m31  -m32  1  0  0\n",
       " -m41     0  0  1  0\n",
       " -m51     0  0  0  1"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E1 * E(3,2,5)  # Why does this have a simple looking answer?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{bmatrix}1&0&0&0&0\\\\m_{21}&1&0&0&0\\\\m_{31}&m_{32}&1&0&0\\\\m_{41}&0&0&1&0\\\\m_{51}&0&0&0&1\\end{bmatrix}"
      ],
      "text/plain": [
       "5×5 Array{SymPy.Sym,2}:\n",
       "   1    0  0  0  0\n",
       " m21    1  0  0  0\n",
       " m31  m32  1  0  0\n",
       " m41    0  0  1  0\n",
       " m51    0  0  0  1"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "inv(E1) * inv(E(3,2,5)) # Why does this have an even slightly simpler looking answer?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{bmatrix}1&0&0&0&0\\\\- m_{21}&1&0&0&0\\\\m_{21} m_{32} - m_{31}&- m_{32}&1&0&0\\\\- m_{41}&0&0&1&0\\\\- m_{51}&0&0&0&1\\end{bmatrix}"
      ],
      "text/plain": [
       "5×5 Array{SymPy.Sym,2}:\n",
       "             1     0  0  0  0\n",
       "          -m21     1  0  0  0\n",
       " m21*m32 - m31  -m32  1  0  0\n",
       "          -m41     0  0  1  0\n",
       "          -m51     0  0  0  1"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    " E(3,2,5) * E1 # Why doesn't this have a simple looking answer?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{bmatrix}1&0&0&0&0\\\\0&1&0&0&0\\\\0&- m_{32}&1&0&0\\\\0&- m_{42}&0&1&0\\\\0&- m_{52}&0&0&1\\end{bmatrix}"
      ],
      "text/plain": [
       "5×5 Array{SymPy.Sym,2}:\n",
       " 1     0  0  0  0\n",
       " 0     1  0  0  0\n",
       " 0  -m32  1  0  0\n",
       " 0  -m42  0  1  0\n",
       " 0  -m52  0  0  1"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E2 = prod( E(i,2,5) for i=3:5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{bmatrix}1&0&0&0&0\\\\0&1&0&0&0\\\\0&0&1&0&0\\\\0&0&- m_{43}&1&0\\\\0&0&- m_{53}&0&1\\end{bmatrix}"
      ],
      "text/plain": [
       "5×5 Array{SymPy.Sym,2}:\n",
       " 1  0     0  0  0\n",
       " 0  1     0  0  0\n",
       " 0  0     1  0  0\n",
       " 0  0  -m43  1  0\n",
       " 0  0  -m53  0  1"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E3 = prod( E(i,3,5) for i=4:5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{bmatrix}1&0&0&0&0\\\\0&1&0&0&0\\\\0&0&1&0&0\\\\0&0&0&1&0\\\\0&0&0&- m_{54}&1\\end{bmatrix}"
      ],
      "text/plain": [
       "5×5 Array{SymPy.Sym,2}:\n",
       " 1  0  0     0  0\n",
       " 0  1  0     0  0\n",
       " 0  0  1     0  0\n",
       " 0  0  0     1  0\n",
       " 0  0  0  -m54  1"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E4 = E(5,4,5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{bmatrix}1&0&0&0&0\\\\- m_{21}&1&0&0&0\\\\- m_{31}&- m_{32}&1&0&0\\\\- m_{41}&- m_{42}&- m_{43}&1&0\\\\- m_{51}&- m_{52}&- m_{53}&- m_{54}&1\\end{bmatrix}"
      ],
      "text/plain": [
       "5×5 Array{SymPy.Sym,2}:\n",
       "    1     0     0     0  0\n",
       " -m21     1     0     0  0\n",
       " -m31  -m32     1     0  0\n",
       " -m41  -m42  -m43     1  0\n",
       " -m51  -m52  -m53  -m54  1"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E1 * E2 * E3 * E4 # Why is this simple?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{bmatrix}1&0&0&0&0\\\\m_{21}&1&0&0&0\\\\m_{31}&m_{32}&1&0&0\\\\m_{41}&m_{42}&m_{43}&1&0\\\\m_{51}&m_{52}&m_{53}&m_{54}&1\\end{bmatrix}"
      ],
      "text/plain": [
       "5×5 Array{SymPy.Sym,2}:\n",
       "   1    0    0    0  0\n",
       " m21    1    0    0  0\n",
       " m31  m32    1    0  0\n",
       " m41  m42  m43    1  0\n",
       " m51  m52  m53  m54  1"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L = inv(E1) * inv(E2)  * inv(E3) * inv(E4) # Why is this simple?  This is the L Matrix!!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{bmatrix}1&0&0&0&0\\\\- m_{21}&1&0&0&0\\\\m_{21} m_{32} - m_{31}&- m_{32}&1&0&0\\\\- m_{21} \\left(m_{32} m_{43} - m_{42}\\right) + m_{31} m_{43} - m_{41}&m_{32} m_{43} - m_{42}&- m_{43}&1&0\\\\- m_{21} \\left(- m_{32} \\left(m_{43} m_{54} - m_{53}\\right) + m_{42} m_{54} - m_{52}\\right) - m_{31} \\left(m_{43} m_{54} - m_{53}\\right) + m_{41} m_{54} - m_{51}&- m_{32} \\left(m_{43} m_{54} - m_{53}\\right) + m_{42} m_{54} - m_{52}&m_{43} m_{54} - m_{53}&- m_{54}&1\\end{bmatrix}"
      ],
      "text/plain": [
       "5×5 Array{SymPy.Sym,2}:\n",
       "                                                                                 1  …              0     0  0\n",
       "                                                                              -m21                 0     0  0\n",
       "                                                                     m21*m32 - m31                 1     0  0\n",
       "                                              -m21*(m32*m43 - m42) + m31*m43 - m41              -m43     1  0\n",
       " -m21*(-m32*(m43*m54 - m53) + m42*m54 - m52) - m31*(m43*m54 - m53) + m41*m54 - m51     m43*m54 - m53  -m54  1"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E4 * E3 * E2 * E1 # Why is this a mess?  Good thing we don't need it in Gaussian Elimination"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{bmatrix}1&0&0&0&0\\\\- m_{21}&1&0&0&0\\\\m_{21} m_{32} - m_{31}&- m_{32}&1&0&0\\\\m_{21} m_{42} - m_{41} - m_{43} \\left(m_{21} m_{32} - m_{31}\\right)&m_{32} m_{43} - m_{42}&- m_{43}&1&0\\\\m_{21} m_{52} - m_{51} - m_{53} \\left(m_{21} m_{32} - m_{31}\\right) - m_{54} \\left(m_{21} m_{42} - m_{41} - m_{43} \\left(m_{21} m_{32} - m_{31}\\right)\\right)&m_{32} m_{53} - m_{52} - m_{54} \\left(m_{32} m_{43} - m_{42}\\right)&m_{43} m_{54} - m_{53}&- m_{54}&1\\end{bmatrix}"
      ],
      "text/plain": [
       "5×5 Array{SymPy.Sym,2}:\n",
       "                                                                               1  …              0     0  0\n",
       "                                                                            -m21                 0     0  0\n",
       "                                                                   m21*m32 - m31                 1     0  0\n",
       "                                             m21*m42 - m41 - m43*(m21*m32 - m31)              -m43     1  0\n",
       " m21*m52 - m51 - m53*(m21*m32 - m31) - m54*(m21*m42 - m41 - m43*(m21*m32 - m31))     m43*m54 - m53  -m54  1"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "inv(L)  # Right this is the inv(L) , not how the inverse of a triangular matrix isn't so pretty in general"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How do we solve Ly=b for the unknown y given b?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{bmatrix}1&0&0&0&0\\\\m_{21}&1&0&0&0\\\\m_{31}&m_{32}&1&0&0\\\\m_{41}&m_{42}&m_{43}&1&0\\\\m_{51}&m_{52}&m_{53}&m_{54}&1\\end{bmatrix}"
      ],
      "text/plain": [
       "5×5 Array{SymPy.Sym,2}:\n",
       "   1    0    0    0  0\n",
       " m21    1    0    0  0\n",
       " m31  m32    1    0  0\n",
       " m41  m42  m43    1  0\n",
       " m51  m52  m53  m54  1"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that the exact answer gets messy:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "b (generic function with 1 method)"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b(i) = symbols(\"b$i\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{bmatrix}b_{1}\\\\b_{2}\\\\b_{3}\\\\b_{4}\\\\b_{5}\\end{bmatrix}"
      ],
      "text/plain": [
       "5-element Array{SymPy.Sym,1}:\n",
       " b1\n",
       " b2\n",
       " b3\n",
       " b4\n",
       " b5"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[b(i) for i=1:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{bmatrix}b_{1}\\\\- b_{1} m_{21} + b_{2}\\\\- b_{1} m_{31} + b_{3} - m_{32} \\left(- b_{1} m_{21} + b_{2}\\right)\\\\- b_{1} m_{41} + b_{4} - m_{42} \\left(- b_{1} m_{21} + b_{2}\\right) - m_{43} \\left(- b_{1} m_{31} + b_{3} - m_{32} \\left(- b_{1} m_{21} + b_{2}\\right)\\right)\\\\- b_{1} m_{51} + b_{5} - m_{52} \\left(- b_{1} m_{21} + b_{2}\\right) - m_{53} \\left(- b_{1} m_{31} + b_{3} - m_{32} \\left(- b_{1} m_{21} + b_{2}\\right)\\right) - m_{54} \\left(- b_{1} m_{41} + b_{4} - m_{42} \\left(- b_{1} m_{21} + b_{2}\\right) - m_{43} \\left(- b_{1} m_{31} + b_{3} - m_{32} \\left(- b_{1} m_{21} + b_{2}\\right)\\right)\\right)\\end{bmatrix}"
      ],
      "text/plain": [
       "5-element Array{SymPy.Sym,1}:\n",
       "                                                                                                                                                              b1\n",
       "                                                                                                                                                    -b1*m21 + b2\n",
       "                                                                                                                               -b1*m31 + b3 - m32*(-b1*m21 + b2)\n",
       "                                                                                     -b1*m41 + b4 - m42*(-b1*m21 + b2) - m43*(-b1*m31 + b3 - m32*(-b1*m21 + b2))\n",
       " -b1*m51 + b5 - m52*(-b1*m21 + b2) - m53*(-b1*m31 + b3 - m32*(-b1*m21 + b2)) - m54*(-b1*m41 + b4 - m42*(-b1*m21 + b2) - m43*(-b1*m31 + b3 - m32*(-b1*m21 + b2)))"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y = L\\[b(i) for i=1:5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Forward Substitution \n",
    "but if you compute the answer the obvious way, it is really straightforward."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$b_{1}$$"
      ],
      "text/plain": [
       "b1"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y[1] = b(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "L[2,1]*y[1] + y[2] = b(2)  implies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$- b_{1} m_{21} + b_{2}$$"
      ],
      "text/plain": [
       "-b1*m21 + b2"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y[2] = b(2) - L[2,1]*y[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "L[3,1]*y[1] + L[3,2]*y[2] + y[3] = b(3)  implies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$- b_{1} m_{31} + b_{3} - m_{32} \\left(- b_{1} m_{21} + b_{2}\\right)$$"
      ],
      "text/plain": [
       "-b1*m31 + b3 - m32*(-b1*m21 + b2)"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y[3] = b(3) - L[3,1]*y[1] - L[3,2]*y[2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Putting this all together"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "y = [b(i) for i=1:5]\n",
    "for i = 1:5, j=1:i-1  \n",
    "    y[i] -= L[i,j]*y[j]\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{bmatrix}b_{1}\\\\- b_{1} m_{21} + b_{2}\\\\- b_{1} m_{31} + b_{3} - m_{32} \\left(- b_{1} m_{21} + b_{2}\\right)\\\\- b_{1} m_{41} + b_{4} - m_{42} \\left(- b_{1} m_{21} + b_{2}\\right) - m_{43} \\left(- b_{1} m_{31} + b_{3} - m_{32} \\left(- b_{1} m_{21} + b_{2}\\right)\\right)\\\\- b_{1} m_{51} + b_{5} - m_{52} \\left(- b_{1} m_{21} + b_{2}\\right) - m_{53} \\left(- b_{1} m_{31} + b_{3} - m_{32} \\left(- b_{1} m_{21} + b_{2}\\right)\\right) - m_{54} \\left(- b_{1} m_{41} + b_{4} - m_{42} \\left(- b_{1} m_{21} + b_{2}\\right) - m_{43} \\left(- b_{1} m_{31} + b_{3} - m_{32} \\left(- b_{1} m_{21} + b_{2}\\right)\\right)\\right)\\end{bmatrix}"
      ],
      "text/plain": [
       "5-element Array{SymPy.Sym,1}:\n",
       "                                                                                                                                                              b1\n",
       "                                                                                                                                                    -b1*m21 + b2\n",
       "                                                                                                                               -b1*m31 + b3 - m32*(-b1*m21 + b2)\n",
       "                                                                                     -b1*m41 + b4 - m42*(-b1*m21 + b2) - m43*(-b1*m31 + b3 - m32*(-b1*m21 + b2))\n",
       " -b1*m51 + b5 - m52*(-b1*m21 + b2) - m53*(-b1*m31 + b3 - m32*(-b1*m21 + b2)) - m54*(-b1*m41 + b4 - m42*(-b1*m21 + b2) - m43*(-b1*m31 + b3 - m32*(-b1*m21 + b2)))"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Back Substitution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For upper triangular matrices with arbitrary diagonal there is an analagous algorithm known as \"back substitution.\"\n",
    "Since the diagonal is not 1, there is one more divide by the diagonal in the obvious alagorithm.\n",
    "In the solution to Ax=b, the pivots are on the diagonal and we divide by those."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.6.0",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
